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An accurate and efficient numerical simulation approach to electromagnetic wave scattering from 
two-dimensional, randomly rough, penetrable surfaces is presented. The use of the Miiller equations 
and an impedance boundary condition for a two-dimensional rough surface yields a pair of coupled 
two-dimensional integral equations for the sources on the surface in terms of which the scattered 
field is expressed through the Franz formulas. By this approach, we calculate the full angular 
intensity distribution of the scattered field that is due to a finite incident beam of p-polarized light. 
We specifically check the energy conservation (unitarity) of our simulations (for the non-absorbing 
case). Only after a detailed numerical treatment of both diagonal and close-to-diagonal matrix 
elements is the unitarity condition found to be well-satisfied for the non-absorbing case (U > 0.995), 
a result that testifies to the accuracy of our approach. 

PACS numbers: 



The scattering of electromagnetic waves from two- 
dimensional randomly rough penetrable surfaces has 
been studied theoretically for more than five decades. 
The methods used in these studies in recent years, where 
attention has been directed toward multiple-scattering 
phenomena, have been either analytical in nature or com- 
putational. Chief among the former methods has been 
the small-amplitude perturbation theory [THS] , while sev- 
eral different computational methods have been used in 
solving the scattering problem. In the earliest calcula- 
tion of this type [4] the system of coupled inhomogeneous 
integral equations for the tangential components of the 
total electric and magnetic fields on the rough surface 
obtained from scattering theory, was converted into a 
system of inhomogeneous matrix equations by the use of 
the method of moments [5], which was then solved by 
Neumann-Liouville iteration [6] . This is a formally exact 
approach but one that is computationally (and memory) 
intensive. 

In subsequent work on this problem approximate solu- 
tions of the exact integral equations have been sought. In 
the sparse-matrix flat-surface iterative approach [TJE] the 
matrix elements for two close points on the surface are 
treated exactly, while those connecting two widely sep- 
arated points are treated approximately, in an iterative 
solution of the matrix equations. 

Soriano and Saillard [9 have combined the sparse- 
matrix flat-surface iterative approach with an impedance 
approximation [10] to calculate co-polarized and cross- 
polarized bistatic scattering coefficients of aluminum ran- 
domly rough surfaces for comparison with results ob- 
tained from perfectly conducting surfaces. 

An approach that combines the fast multipole 
method [IT] and the sparse-matrix flat-surface iterative 
approach has been developed by Jandhyala et. al [T2] . 



Despite these advances, the calculation of the scatter- 
ing of electromagnetic waves from two-dimensional, pene- 
trable, randomly rough surfaces, remains a computation- 
ally intensive problem, in need of further improvements 
in the methods used to solve it. 

In this paper we use the Franz formulas of electromag- 
netic scattering theory [13] [14] to obtain expressions for 
the amplitude of the electromagnetic field scattered from 
a two-dimensional, rough, metallic or dielectric surface 
in terms of the tangential components of the total elec- 
tric and magnetic fields on the surface. The independent 
elements of these tangential field components satisfy a 
system of four coupled inhomogeneous two-dimensional 
integral equations — the Miiller equations [15] [16] — de- 
rived from Franz formulas. This system of four integral 
equations is reduced to a system of two integral equa- 
tions by the use of an impedance boundary condition at 
a two-dimensional rough metallic surface [17] , and its so- 
lution is used to calculate the mean differential reflection 
coefficient. 

The approach to the scattering of an electromagnetic 
field from a rough metallic or dielectric surface outlined 
here is similar to the approach of Soriano and Saillard [5 
in its use of an impedance boundary condition to reduce 
the number of coupled integral equations that have to be 
solved from four to two. However, there are still impor- 
tant differences between these two approaches. The first 
is that we do not use the sparse-matrix flat-surface iter- 
ative approach: the matrix elements connecting any two 
points on the surface are calculated accurately. More- 
over, those connecting two nearby points are calculated 
more accurately than in the work of Soriano and Saillard. 
The second is that we calculate the full angular inten- 
sity distribution of the scattered field, which allows us 
to check the unitarity (energy conservation) of the scat- 
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FIG. 1: (Color online) The mean differential reflection coef- 
ficients [20] , (dRvp/dQs) (p — > v) for a p-polarized incident 
beam, whose polar angle of incidence is Oq = 20°, as functions 
of the polar scattering angle S for the (a) in-plane and (b) 
out-of-plane scattering. See text for additional parameters. 

tered field (for the non-absorbing case). This enabled the 
identification of critical aspects of the numerical imple- 
mentation that, if not handled properly, could lead to 
erroneous results and/or a significant drop in accuracy. 
This important point seems to have been overlooked in 
previous publications. The third is that although the 
occurrence of hyper-singular kernels is avoided in both 
approaches by the use of the Miiller equations [I5j [16] , 
some differences are found between our resulting matrix 
elements and those of Soriano and Saillard that appear to 
affect the unitarity of the scattered field[24]. Moreover, 
contrary to what was reported in Ref. [15 , we find that 
matrix element terms containing the Green's function of 
the metal also have to be taken into account for some off- 
diagonal elements in order to produce accurate results. 
The fourth is that we do not use the beam decomposi- 
tion method [19] for the incident beam in which a wide 
beam is represented by a weighted sum of shifted narrow 
beams. Instead we use a single wide incident beam. 

The physical system that we consider consists of 
vacuum [e = 1] in the region x 3 > C( x ||) [where 
X|| = (xi,X2,0)] and a non-magnetic metal in the re- 
gion xs < C( x ||) that is characterized by a complex, 
frequency-dependent, dielectric function, e(uj), for which 
Ree(uj) < and Ime(uo) > 0. The surface profile func- 
tion, xs = C( x ||)? i s assumed to constitute a zero-mean, 
Gaussian random process that is a single-valued func- 
tion of X|| and is differentiate with respect to x\ and 
X2 at least twice. It is defined by (C( x ||)) = an d 
C(x||)C(xj|)^ = (5 2 W(|x|| - xj||), where S is the root- 
mean-square roughness, W(-) denotes the (normalized) 
correlation function, and the angle brackets denote an 
average over the ensemble of realizations of C( x ||)- I n this 
work we will use an isotropic Gaussian correlation func- 
tion W(x) = exp(— x 2 / a 2 ) with a the correlation length. 



The starting point of our approach is the Franz formu- 
las of electromagnetic theory (or the dyadic form of Huy- 
gens principle) [131 [14] . By applying them to the vacuum 
region above the metal surface, and letting the observa- 
tion point, x, approach the surface xs = C( x ||)? two cou- 
pled inhomogeneous integral equations for the tangential 
components of the electric and magnetic fields, J#(x||) = 

& x E ( x )U=C(x,|) and J ^( x ll) = & X H ( x )U=C(x„)> re " 
spectively, are obtained, where n denotes the unit vec- 
tor normal to the surface and directed into the vac- 
uum. These equations contain double derivatives of the 
scalar Green's function go( x ? x ') = ex.p[ikoR]/47rR [ko = 
y / 5o(cJ) cj/c, R = |x — x'|] resulting in non-integrable 
hyper- singular kernels that are sources of computational 
difficulties [9] [15j [16] . One way to obtain integrable ker- 
nels, is to combine in a suitable way the two sets of Franz 
formulas obtained separately for the vacuum and metal 
regions so that the resulting integral equations do not 
contain any hyper-singular terms. The resulting equa- 
tions are known as the Miiller integral equations [15], 
and the one satisfied by Jif(xn \ w) reads 

J#(X|||CJ) = J#(X|||cj) inc 

+ pj d 2 x\ [n(x||) x {V x [S(x|x') J H (^)}}] (1) 
- ^ J d 2 x\ [n(x||) x {V x V x [0(x|x') J^xj, |^)]}], 

where the equation satisfied by J#(xj||u;) can be ob- 
tained from duality [14 . In writing Eq. we have 
introduced <7(x|||xj,) = # ( x ||| x j|) - #( x ||l x j|) — the dif- 
ference between the scalar Green's functions for the 
vacuum (subscript 0) and the metal; [A(x|x')] = 
=C(x||);4=C(x(i)' ^ denotes the Cauchy princi- 
ple value of an integral; and Jh( x || \^)inc is defined sim- 
ilarly to J#(x|||co>) but for the incident field. Initially 
the kernel of Eq. ([I]) seems to be hyper-singular. How- 
ever, because the leading term (when R — >• 0) of the 
second derivative of the scalar Green's function is in- 
dependent of medium parameters, the most divergent 
terms of the kernel cancel, rendering it integrable. By 
adopting the impedance boundary condition that re- 
lates the surface currents J#(x|||cj) and J j h-(xn \uj) via 
the (local) impedance tensor (K) [17]: J^(xn|a;)i = 
Kij(x\\ |o;)Jij(x|| [i,j = 1,2], the dependence on 
J#(x|| \uj) can be removed from Eq. |l]). Moreover, the re- 
sulting equation can be converted into a matrix equation 
for the two independent electric surface current compo- 
nents, say, Jjf(xi| \uj)i [i — 1, 2], by the use of the method 
of moments [5] . The resulting linear system is then solved 
by the BiCGStab method [18 and the solution used to 
calculate the mean differential reflection coefficient that 
is averaged over an ensemble of realizations of the surface 
profile function (see Ref. [20] for details). 

On the basis of the integral equation 0, and with 
the use of the impedance boundary condition, we have 
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performed numerical simulations for a p-polarized inci- 
dent beam of wavelength A = 0.6328 fim that is scat- 
tered from a Gaussian randomly rough silver surface. At 
this wavelength e(u) = -16.00 + zl.088 [21]. The sur- 
face was characterized by a root-mean-square roughness 
of 5 = A/4 and a correlation length a = A/2. In the sim- 
ulations it was assumed to cover an area of 16A x 16A, 
and the discretization interval was A = A/7 for both the 
x\- and ^-directions. 

Figure [I] presents the mean differential reflection co- 
efficients as functions of the polar scattering angle S 
for the in-plane [Fig. [TJa)] and out-of-plane (cj) s = 
±90°) [Fig. [TJb)] scattered light due to a p-polarized 
Gaussian beam of width w = 4A incident on the sur- 
face at a polar angle 0o = 20°. For the same parameters, 
Figs. [2] depict the full angular intensity distribution of 
the incident p-polarized light that is scattered into p- and 
5-polarized light (polarization not recorded) [Fig. 2 
p-polarization [Fig. j2^b)] ; and s-polarization [Fig. 2 
The number of surface realizations used to obtain these 
results was = 5000. The simulations required for ev- 
ery surface realization 96 CPU seconds (on a 2.67GHz 
Intel i7 CPU) for each angle of incidence when calculat- 
ing the scattered field on a 100 x 100 grid. The peaks 
observed in Figs, [I] and [2] at angular positions S = — 6$ 
(and <j) s = (/>o+7r) are due to the enhanced backscattering 
phenomenon, a multiple scattering effect [22 . The en- 
ergy fraction of the incident light that is scattered by the 
surface was 94.7%, compared to 96.9% as predicted from 
the Fresnel coefficient of the corresponding flat interface 
scattering system. All the light scattered by the surface 
was essentially incoherent (diffuse) (about 99.98%). 

In order to evaluate the accuracy of the simulations 
and to perform a self-consistency check of our approach, 
we have performed simulations using the parameters 
given above with the exception that the metal was as- 
sumed to be non-absorbing, i.e. we artificially put 
lme(u) = 0. Under this assumption, the total time- 
averaged power fluxes of the incident and scattered fields 
have to be equal, or in other word, one should require en- 
ergy conservation (or equivalently unitarity of the scat- 
tered field, U = 1 where U denotes the fraction of the 
incident power flux that is scatted by the rough metal 
surface. It should be stressed that energy conservation 
is only a necessary, but not sufficient condition to guar- 
antee the correctness of the simulation results for non- 
absorbing systems. It is still, however, a rather useful 
and non-trivial condition that can assist in detecting in- 
accuracies of the simulation approach as well as potential 
implementation errors. For the parameters used in the 
simulations reported in this work, we found U > 0.995 
for "non-absorbing" silver [s(uj) = —16.00], a result that 
testifies to the accuracy of our approach. 

In order to achieve such good unitarity values, it 
turned out that great care had to be exercised when han- 
dling the matrix elements of the integral equation kernel 
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FIG. 2: (Color online) The same as Fig. [I] but now showing 
the full angular intensity distribution of the scattered field, 
(a) p — »• polarization not recorded; (b) p — »• p; and (c) p — >> s. 



that were on, or close to, the diagonal. Soriano and Sail- 
lard j9] have previously pointed out one way of handling 
the diagonal matrix elements that contain the singular- 
ity (at X|| = xj|) of the Green's function by separating 
it into two parts: one for which the integrand is singu- 
lar but integrable and is handled analytically, and an- 
other that is regular and is handled numerically. These 
authors were not able within their approach to check 
the energy conservation of their calculations. We have 
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found, however, that in order to achieve good results 
for the energy conservation, also close-to diagonal ma- 
trix elements (in addition to those on the diagonal) have 
to be treated accurately even if the off-diagonal matrix 
elements are regular everywhere. These findings some- 
what resemble results reported for volume integral equa- 
tions where also close to diagonal volume-elements had 
to be handled with higher accuracy then more distant 
matrix elements [23]. For instance the use of the mid- 
point method for evaluating all off-diagonal matrix ele- 
ments (as in Ref. [9 j) and a more accurate method for 
the diagonal elements, would, for the surface parameters 
assumed here, result in about 16.4% (U = 0.834) of the 
incident energy not being accounted for, a result that 
was found to be more-or-less independent of how accu- 
rately one treated the diagonal elements, or if the surface 
was rough or flat. Moreover, if in addition to the diago- 
nal matrix elements, also the nearest-neighbor elements 
were treated accurately, the amount of energy that was 
not accounted for dropped to 4.9% (U = 0.951). If, addi- 
tionally, also next-nearest neighbor matrix elements were 
treated accurately, the unitarity condition started to be- 
come well satisfied (U > 0.995), and treating accurately 
matrix elements that were even further apart contributed 
only insignificantly to the improvement of the unitarity 
condition [24 . It should be stressed that without per- 
forming the self-consistency check based on energy con- 
servation, which requires the full angular distribution of 
the scattered light being available to us, it has probably 
not been realized that failing to treat close-to-diagonal 
matrix elements accurately could cause inaccuracies in 
the range of 10-20% even for flat interfaces. This is one 
of the main results of this Letter. 

In conclusion, we have presented an accurate and high- 
performance simulation approach for the scattering of 
electromagnetic waves from two-dimensional penetrable 
metallic surfaces based on surface integral techniques. 
By this approach, the scattering of a p-polarized finite 
beam by a two-dimensional, randomly rough, silver sur- 
face was studied in the optical regime, and it gave rise 
to the well-known enhanced backscattering phenomenon. 
Due to the calculation of the full angular intensity dis- 
tribution of the scattered light, it was possible for us 
to evaluate the accuracy of the simulation approach. It 
was found that high-quality simulation results required 
an accurate treatment of both the diagonal and close- 
to-diagonal matrix elements. This latter point seems to 
have been overlooked in previous studies. In this way, we 
were able to obtain results that respect energy conserva- 
tion (unitarity) for the equivalent non-absorbing system, 
something that testifies to the accuracy of our approach. 

The simulation approach presented in this Letter opens 
the door to a direct and detailed comparison of the 
full angle-resolved intensity distributions of the scattered 
light obtained experimentally and theoretically. Addi- 



tionally, the approach provides the tools needed to pre- 
dict the effect of surface roughness on the electromagnetic 
field in the near and far zone of the surface, and also 
to tailor surface structures towards specific applications 
(engineered surfaces). Such issues are of importance in 
numerous applications, like for instance, in the photo- 
voltaic industry where surface roughness in solar cells is 
known to increase the efficiency of the cell, but the op- 
timal surface structure, and the mechanism responsible 
for the increased efficiency, are still unknown [25] . 
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(Smaforsk grant), and NTNU. 
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